require(ggplot2)
require(reshape2)
require(dplyr)
require(tidyr)
require(ncdf4)
require(raster)
require(grid)
require(lubridate)
library(suncalc) 
require(stats)

setwd("~/Desktop/grass_exp/cal/")
source("~/Desktop/tsj_plotting.R")


#spectroradiometer:-----------------------------------------------------------------------
in_dir <- "spectroradiometer_data/2020_Oct_08/"
in_files <- list.files(in_dir)
header_df <- data.frame()
for(in_file in in_files){
  if( strsplit(in_file,"\\.")[[1]][2] == "sed"){
    print(in_file)
    line = ""
    header_df_row <- data.frame(in_file=in_file)
    header_names <- c("in_file")
    con = file(paste0(in_dir,in_file),"r")
    while( line != "Data:"){
      line = readLines(con, n = 1)
      var = strsplit(line,": ")[[1]][1]
      val = strsplit(line,": ")[[1]][2]
      header_df_row <- cbind( header_df_row, data.frame( val, stringsAsFactors = F) )
      header_names <- c( header_names, var)
    }
    close(con)
    names(header_df_row) <- header_names
    header_df_row$Date <- strsplit(header_df_row$Date[1],",")[[1]][1]
    header_df_row$Time <- strsplit(header_df_row$Time[1],",")[[1]][2]
    header_df <- rbind( header_df, header_df_row)
  }
}

header_df$tstamp <- as.POSIXct( paste0( header_df$Date, header_df$Time),
                                format = "%m/%d/%Y%H:%M:%S")

t0 <- as.POSIXct("2020-10-08 00:00:00",tz='US/Eastern')
header_df$footLamberts <- NA

header_df$footLamberts[ (header_df$tstamp > t0 + hours(16) + minutes(20)) &
                          (header_df$tstamp < t0 + hours(16) + minutes(21))] <- 0.5012e3

header_df$footLamberts[ (header_df$tstamp > t0 + hours(16) + minutes(25)) &
                          (header_df$tstamp < t0 + hours(16) + minutes(26))] <- 1.0007e3

header_df$footLamberts[ (header_df$tstamp > t0 + hours(16) + minutes(32)) &
                          (header_df$tstamp < t0 + hours(16) + minutes(33))] <- 1.0011e2

header_df$footLamberts[ (header_df$tstamp > t0 + hours(16) + minutes(36)) &
                          (header_df$tstamp < t0 + hours(16) + minutes(37))] <- 0.5000e2

header_df$footLamberts[ (header_df$tstamp > t0 + hours(16) + minutes(40)) &
                          (header_df$tstamp < t0 + hours(16) + minutes(41))] <- 0

srad_df <- data.frame()
for( i in 1:dim(header_df)[1]){
  if( !is.na(header_df$footLamberts[i])){
    srad_df_i <- read.table(paste0(in_dir,header_df$in_file[i]),
                            skip = length(header_df), 
                            stringsAsFactors = F)
    names(srad_df_i) <- c("wl","rad_ref","rad")
    srad_df_i$footLamberts <- header_df$footLamberts[i]
    srad_df <- rbind( srad_df, srad_df_i)
  }
}
#srad_df <- aggregate(rad~wl+footLamberts,srad_df,"mean")
wls <- unique(srad_df$wl)
rad_cal_df <- data.frame()
for(wl in wls){
  sdf <- srad_df[ srad_df$wl == wl, ]
  foo <- summary(lm(rad~footLamberts,sdf))
  rad_cal_df <- rbind( rad_cal_df, data.frame(wl=wl,
                                              b=foo$coefficients[1],
                                              b_err=foo$coefficients[1,2],
                                              a=foo$coefficients[2],
                                              a_err=foo$coefficients[2,2]))
}

#just the fld wavelengths:
out_band_ideal_wl = 757.5
in_band_ideal_wl = 760.5
wl_out_idx <- which.min( abs( wls - out_band_ideal_wl ) )
wl_in_idx  <- which.min( abs( wls -  in_band_ideal_wl ) )
sdf    <- srad_df[ srad_df$wl %in% c( wls[wl_in_idx], wls[wl_out_idx]), ]
sdf$wl <- as.character(sdf$wl)
ggplot( sdf, aes(x=footLamberts,y=rad,color=wl)) +
  geom_point() +
  theme_bw() 


ggplot(rad_cal_df,aes(x=wl,y=a*1e3))+
  geom_point() +
  geom_errorbar(aes(ymin=(a-a_err)*1e3,ymax=(a+a_err)*1e3)) +
  xlim(730,790) +
  theme_tsj() +
  labs(x="Wavelength [nm]",
       y=bquote("[mW"~m^{-2}~nm^{-1}~sr^{-1}~"] / fL" )) +
  ylim(.1,.12)

ggplot(rad_cal_df,aes(x=wl,y=b*1e3))+
  geom_point() +
  geom_errorbar(aes(ymin=(b-b_err)*1e3,ymax=(b+b_err)*1e3)) +
  xlim(730,790) +
  theme_tsj() +
  labs(x="Wavelength [nm]",
       y=bquote("[mW"~m^{-2}~nm^{-1}~sr^{-1}~"]" ))

#Cameron Data:------------------------------------------------------------------

qe <- "QEP02110"

dates         <- c("20210720")
data_dir      <- "rad_cal_spectra/"
dark_cal_file <- "QEP02723_dark_calibration_raw_values_20200428-1416.csv"
foot_lamberts <- c(        1.0115e3,         0.797e3,         0.571e3,         1.401e3,         1.656e3,         2.491e3)
start_times   <- c("20210720 15:45","20210720 15:50","20210720 15:56","20210720 16:10","20210720 16:18","20210720 16:26")
end_times     <- c("20210720 15:47","20210720 15:52","20210720 15:58","20210720 16:12","20210720 16:20","20210720 16:27")
modes         <- c(        "normal",        "normal",        "normal",        "normal",        "normal",        "normal")
#----------------------------------------------------------------------------

stimes <- as.POSIXct( start_times, format = "%Y%m%d %H:%M", tz="UTC") #+ 4*60*60
etimes <- as.POSIXct( end_times,   format = "%Y%m%d %H:%M", tz="UTC") #+ 4*60*60

thermal_df <- data.frame()
df <- data.frame()
spec <- NULL
inttime <- NULL
dtime <- NULL
for( date in dates){
  cam_data_file <- paste0(data_dir, "cameron_",qe,"_",date,"_sphere.nc" )
  nc <- nc_open(cam_data_file)
  
  thermal_df <- rbind( thermal_df, 
                       data.frame( date = date,
                                   time = ncvar_get(nc,"thermal/time"),
                                   spec_temp = ncvar_get(nc, "thermal/spectrometer_temperature" ),
                                   board_temp = ncvar_get(nc, "thermal/board_temperature" ),
                                   tec_power = ncvar_get(nc, "thermal/TEC_power" ) ) ) 

  target    <- "fixed"
  wl        <- ncvar_get(nc,paste0(target,"/wavelength") )
  spec      <- cbind( spec,  ncvar_get( nc,paste0(target,"/spectra") ) )
  dtime     <- c( dtime, as.POSIXct("1970-01-01", tz="UTC") + ncvar_get(nc,paste0(target,"/time")))
  inttime   <- c( inttime, ncvar_get(nc,paste0(target,"/integration_time")))
}

dtime <- as.POSIXct( dtime, origin = "1970-01-01")
timing_df <- data.frame( dtime = dtime, inttime=inttime )

ggplot( timing_df, aes(x=dtime,y=inttime)) +
  geom_line() #+
  #xlim(stimes[1],etimes[length(etimes)])
  

timing_df$mode <- NaN
timing_df$fL <- NaN
for( i in 1:length(foot_lamberts)){
  timing_df$fL[ (timing_df$dtime > stimes[i]) & (timing_df$dtime < etimes[i]) ] <- foot_lamberts[i]
  timing_df$mode[ (timing_df$dtime > stimes[i]) & (timing_df$dtime < etimes[i]) ] <- modes[i]
}
spec <- spec[,!is.nan(timing_df$fL)]
inttime <- inttime[ !is.nan(timing_df$fL)]
fL <- timing_df$fL[ !is.nan(timing_df$fL)]
mode <- timing_df$mode[ !is.nan(timing_df$fL)]

ggplot( timing_df, aes(x=dtime,y=inttime,color=as.character(fL))) +
  geom_point() +
  theme_tsj() +
  ylim(0,2e6)
  

#read in dark calibration:------------------------------
dark_df <- read.csv(dark_cal_file)
fit_in  <- summary ( lm(  DN_in ~ int_time, data=dark_df) )
#fit_out <- lm( DN_out ~ int_time, data=dark_df)

b <- fit_in$coefficients[1]
a <- fit_in$coefficients[2]

spec2 <- spec - (a*inttime + b)
spec3 <- spec2 %*% diag(1/inttime)

rad_a <- approx( rad_cal_df$wl, rad_cal_df$a, wl )$y
rad_b <- approx( rad_cal_df$wl, rad_cal_df$b, wl )$y

true_rad <- matrix(rad_a)%*%t(matrix(fL)) + rad_b

sp_df <- data.frame()
foo <- data.frame()

for( mode_oi in unique(mode)){
  print(mode_oi)
  for(i in 1:length(wl)){
    fit <- summary( lm( spec3[i,mode==mode_oi]~true_rad[i,mode==mode_oi]) )
    sp_df <- rbind( sp_df, data.frame(mode=mode_oi,
                                      wl=wl[i],
                                      a=fit$coefficients[2],
                                      b=fit$coefficients[1],
                                      a_err=fit$coefficients[2,2],
                                      b_err=fit$coefficients[1,2]))
    
    foo <- rbind( foo, data.frame(mode=mode_oi,
                                      wl=wl[i],
                                      cal=true_rad[i,mode==mode_oi]/spec3[i,mode==mode_oi]))
  }
}

#just the fld wavelengths:
library(ggpubr)

out_band_ideal_wl = 757.5
in_band_ideal_wl = 760.5
wl_out_idx <- which.min( abs( wl - out_band_ideal_wl ) )
wl_in_idx  <- which.min( abs( wl -  in_band_ideal_wl ) )
sdf <- data.frame( DN = spec3[wl_out_idx,], rad = true_rad[wl_out_idx,], wl=wl[wl_out_idx])
sdf <- rbind( sdf, data.frame( DN = spec3[wl_in_idx,], rad = true_rad[wl_in_idx,], wl=wl[wl_in_idx]) )
sdf$wl <- round(sdf$wl,2)
sdf$wl <- as.character(sdf$wl)
ggplot( sdf, aes(x=DN,y=rad,color=wl)) +
  geom_point() +
  theme_bw() +
  labs(x = bquote("Corrected D.N.  ["~mu*s^{-1}~"]"),
       y = bquote("Radiance  [ W"~m^{-2}~nm^{-1}~sr^{-1}~"]"),
       color="Wavelength") +
  geom_smooth(method = "lm", se = FALSE, alpha=.5) +
  stat_regline_equation( aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")) )


sdf <- data.frame( int_time = dark_df$int_time, DN = dark_df$DN_in, wl=wl[wl_in_idx])
sdf <- rbind( sdf, data.frame( int_time = dark_df$int_time, DN = dark_df$DN_out, wl=wl[wl_out_idx]))
sdf$wl <- round(sdf$wl,2)
sdf$wl <- as.character(sdf$wl)
ggplot(sdf,aes(x=int_time,y=DN,color=wl)) +
  geom_point(alpha = .5) +
  theme_bw() +
  geom_smooth(method = "lm", se = FALSE, alpha=.5) +
  stat_regline_equation( aes(label =  paste(..eq.label.., ..adj.rr.label.., sep = "~~~~")) ) +
  labs(x="Integration Time [ microseconds ]", 
       y="D.N.")



ggplot(sp_df, aes(x=wl,y=a)) +
  geom_line() +
  #xlim(730,785) +
  #ylim(10,28) +
  theme_tsj() +
  labs(x="Wavelength [nm]",y = "rad / ( DN / s)" )

ggplot(sp_df, aes(x=wl,y=b)) +
  geom_point() +
  xlim(730,785) 
  #ylim(0,30)

write.csv(sp_df,paste0("Rad_Calibration_",qe,"_",dates[1],".csv"))

foo2 <- foo[ foo$mode != "old_fiber", ]
foo2 <- foo2[ foo2$mode != "new_fiber5", ]
foo2 <- foo2[ foo2$mode != "new_fiber6", ]

foo2 <- foo2[ foo2$wl ==  unique(foo$wl[50]) , ]


foo2$mode[ foo2$mode =="old_fiber2"] <- "old_fiber"



ggplot( foo, aes(x=mode,y=cal,color=mode)) +
  geom_boxplot() +
  #geom_line(aes(y=rad, color="rad")) +
  #geom_point(aes(y=rad, color="rad")) +
  #xlim(731,780) +
  #ylim(.7,1.1) +
  theme_tsj() +
  labs( x= "Wavelength [nm]", y= "rad / ( DN / s)" )



foo3 <- aggregate(cal~wl+mode,foo2,"mean")

foo3$rel <- foo3$cal/foo3$cal[ foo3$mode == "new_fiber"]

foo3 <- foo3[ foo3$mode != "new_fiber3", ]

ggplot(foo3,aes(x=wl,y=100*rel,color=mode)) +
  geom_point() +
  xlim(745,758) +
  theme_tsj() +
  ylim(85,130) +
  labs(x="Wavelength [nm]",y="Relative Calibration [%]")

for(m in unique(foo3$mode)){
  print(m)
  print( summary( lm(rel~wl,foo3[foo3$mode==m,]) ) )
}
#write.csv(sp_df,paste0("Rad_Calibration_",qe,"_",dates[1],".csv"))


ggplot(sub_df,aes(x=wl,y=a))+
  geom_line() +
  xlim(651,860) +
  ylim(10,28) +
  theme_tsj() +
  #facet_wrap(~sensor,scales="free_y",nrow=2) +
  labs(x="Wavelength [nm]",y=bquote("Radiance / (DN/"~mu~"s)"))

i <- 100
sp_df <- data.frame( wl=wl, raw=spec[,i],dk=spec2[,i],dkit=spec3[,i],rad=true_rad[,i])

sp_df$cal <- sp_df$rad/sp_df$dkit

ggplot( sp_df, aes(x=wl)) +
  geom_line(aes(y=dkit,color="DN/us")) +
  #geom_line(aes(y=rad, color="rad")) +
  #xlim(731,780) +
  #xlim(745,758) +
  ylim(.75,1) +
  theme_tsj()

foo <- sp_df[ sp_df$wl ==  unique(sp_df$wl[50]) , ]

ggplot( foo, aes(x=mode,y=mu)) +
    geom_boxplot() +
    #geom_line(aes(y=rad, color="rad")) +
    #geom_point(aes(y=rad, color="rad")) +
    #xlim(731,780) +
    ylim(1,1.2) +
    theme_tsj() +
    labs( x= "Wavelength [nm]", y= "rad / ( DN / s)" )
  

#-------------------------------




ggplot(df, aes(x=time,y=inttime)) +
  geom_point()

df$footLamberts <- NA
for(i in 1:length(foot_lamberts)){
  df$footLamberts[ (df$time > stimes[i]) & (df$time < etimes[i])] <- foot_lamberts[i]
}

df <- df[ !is.na(df$footLamberts), ]
df <- df[ df$footLamberts != 0, ]

#df$sample_no <- 0
#for( target in targets){
#  sdf2 <- df[ df$target == target , ]
#  df$sample_no[ df$target == target] <- 1 + c(0, cumsum( diff(sdf2$time) > 3 ) )
#}

#noise_df <- aggregate( .~ sample_no+target+date,df,"sd" )
#sdf2 <- aggregate( .~ sample_no+target+date,df,"median" )
#n_df <- aggregate(.~ sample_no+target+date,df,"length")
#sdf2$raw_i_sd <- noise_df$raw_i
#sdf2$n_samples <- n_df$raw_i




ggplot(df) +
  geom_histogram(aes(x=raw_i, fill="in_well"),alpha=.5) +
  geom_histogram(aes(x=raw_o, fill="out_well"),alpha=.5) +
  facet_wrap(~footLamberts,ncol=1,scales="free_y") 

ggplot(df,aes(x=time,y=inttime/1e6)) +
  geom_line() +
  facet_wrap(~footLamberts,scales="free")



thermal_df$time <- as.POSIXct("1970-01-01", tz="UTC") + thermal_df$time
thermal_df$hour <- hour(thermal_df$time) + minute(thermal_df$time)/60.0 + second(thermal_df$time)/3600.0


#read in dark calibration:------------------------------
dark_df <- read.csv(dark_cal_file)
fit_in  <- lm(  DN_in ~ int_time, data=dark_df)
fit_out <- lm( DN_out ~ int_time, data=dark_df)
 
df$DK_i <- df$raw_i - ( df$inttime*fit_in$coefficients[2]  + fit_in$coefficients[1] )
df$DK_o <- df$raw_o - ( df$inttime*fit_out$coefficients[2] + fit_out$coefficients[1] )

#DN/ms:
df$i  <- df$DK_i/(df$inttime*1e-3)
df$o <- df$DK_o/(df$inttime*1e-3)



#df$sample_no <- 0
#for( target in targets){
#  sdf2 <- df[ df$target == target , ]
#  df$sample_no[ df$target == target] <- 1 + c(0, cumsum( diff(sdf2$time) > 3 ) )
#}
#noise_df <- aggregate( .~ sample_no+target+date,df,"sd" )
#sdf2 <- aggregate( .~ sample_no+target+date,df,"median" )
#n_df <- aggregate(.~ sample_no+target+date,df,"length")
#sdf2$i_sd <- noise_df$i
#sdf2$o_sd <- noise_df$o
#sdf2$n_samples <- n_df$i

ggplot(df) +
  geom_histogram(aes(x=i, fill="in_well"),alpha=.5) +
  geom_histogram(aes(x=o, fill="out_well"),alpha=.5) +
  facet_wrap(~footLamberts,ncol=1,scales="free") 

ggplot(df) +
  geom_histogram(aes(x=i, fill="in well"),alpha=.5) +
  geom_histogram(aes(x=o, fill="out well"),alpha=.5) +
  facet_wrap(~footLamberts,ncol=2,scales="free") +
  theme_tsj() +
  labs(x="Dark Corrected DN/ms",y="count",fill="") +
  theme(legend.position = "bottom")

ggplot(df,aes(x=footLamberts)) +
  theme_tsj() +
  geom_point(aes(y=i,color="in_well")) +
  geom_point(aes(y=o,color="out_well")) +
  geom_smooth(aes(y=i,color="in_well"),method = lm,size=.5,se = F) +
  geom_smooth(aes(y=o,color="out_well"),method = lm,size=.5,se=F) +
  labs(x="Sphere Output [foot-lamberts]",
       y="Dark Corrected DN/ms",
       title="Sphere vs. Cameron",
       color="")

ggplot(df,aes(x=footLamberts)) +
  theme_tsj() +
  geom_boxplot(aes(y=inttime*1e-6,group=footLamberts)) +
  labs(x="Sphere Output [foot-lamberts]",
       y="Integration Time [sec]",
       title="Integration Times")


#if sphere is calibrated:
df$real_i <-  in_well_sphere_a*df$footLamberts + in_well_sphere_b
df$real_o <- out_well_sphere_a*df$footLamberts + out_well_sphere_b

fit <- lm( i ~ real_i, df)
in_well_b <- fit$coefficients[1]
in_well_a <- fit$coefficients[2]  

fit <- lm( o ~ real_o, df)
out_well_b <- fit$coefficients[1]
out_well_a <- fit$coefficients[2]

ggplot(df) +
  geom_point(aes(y=i,x=real_i*1e3,color="in well")) +
  geom_point(aes(y=o,x=real_o*1e3,color="out well")) +
  geom_smooth(aes(y=i,x=real_i*1e3,color="in well"),method = lm,size=.5,se=F) +
  geom_smooth(aes(y=o,x=real_o*1e3,color="out well"),method = lm,size=.5,se=F) +
  theme_tsj() +
  labs(x = bquote("Dark Corrected Counts [DN"~s^{-1}~"]"),
       y = bquote("Radiance [W"~m^{-2}~nm^{-1}~sr^{-1}~"]"),
       color="") +
  geom_text(aes(x=50,
                y=10,
                label = paste0(
                  "y = ",
                  as.character(round( in_well_a*1e-3, 4)),
                  "x + ",
                  as.character(round( in_well_b*1e-3, 3))),
                color="in well")) +
  geom_text(aes(x=50,
                y=20,
                label = paste0(
                  "y = ",
                  as.character(round( out_well_a*1e-3, 4)),
                  "x + ",
                  as.character(round( out_well_b*1e-3, 3))),
                color="out well"))

ggplot(df) +
  geom_point(aes(y=i,x=real_i*1e3) ) +
  geom_smooth(aes(y=i,x=real_i*1e3,color="in well"),method = lm,size=.5,se=F) +


plt_thermal_tec <- ggplot( thermal_df, aes(x=hour,y=tec_power)) + 
  geom_line() +
  theme_tsj() + 
  ylim(-100,100) +
  labs(x="Hour of Day [UTC]",y="TEC Heating [Percent]") +
  facet_wrap(~date) +
  xlim(12,24)

plt_thermal_temp <- ggplot( thermal_df, aes(x=time)) +
  geom_point(aes(y=spec_temp,color="Spectrometer Bay")) +
  geom_line(aes(y=board_temp,color="TEC Controller")) +
  theme_tsj() +
  labs(x="Time [UTC]", y="Temperature [ Deg C ]", color="") +
  xlim( as.POSIXct("2020-10-08 14:00:00",tz="UTC"), 
        as.POSIXct("2020-10-08 15:30:00",tz="UTC") ) +
  theme(legend.position = "bottom")

plt_int_time <- ggplot( df, aes(x=hour, y=inttime*1e-6, color=target)) +
  geom_line(alpha=.8) +
  theme_tsj() + 
  labs(x= "Time [UTC]", y= "Int. Time [sec]",title=paste0("Integration Times: ")) +
  #ylim(0,5) +
  facet_wrap(~date,ncol=1)




header_df <- header_df[ !is.na(header_df$footLamberts), ]

ggplot( header_df, aes(x=footLamberts)) +
  theme_tsj() +
  geom_point(aes(y=in_well,color="in well")) +
  geom_point(aes(y=out_well,color="out well")) +
  geom_smooth(aes(y=in_well,color="in well"),method = lm,size=.5,se=F) +
  geom_smooth(aes(y=out_well,color="out well"),method = lm,size=.5,se=F) +
  labs(x = "Sphere Output [foot-lamberts]",
       y = bquote("Radiance [W"~m^{-2}~nm^{-1}~sr^{-1}~"]"),
       color="")

ggplot(header_df) +
  geom_histogram(aes(x=in_well, fill="in_well"),alpha=.5) +
  geom_histogram(aes(x=out_well, fill="out_well"),alpha=.5) +
  facet_wrap(~footLamberts,ncol=2,scales="free") 

fit <- lm(in_well~footLamberts,header_df)
in_well_b <- fit$coefficients[1]
in_well_a <- fit$coefficients[2]

fit <- lm(out_well~footLamberts,header_df)
out_well_b <- fit$coefficients[1]
out_well_a <- fit$coefficients[2]

df$real_i <-  in_well_a*df$footLamberts + in_well_b
df$real_o <- out_well_a*df$footLamberts + out_well_b

ggplot(df) +
  geom_point(aes(x=i,y=real_i*1e3,color="in well")) +
  geom_point(aes(x=o,y=real_o*1e3,color="out well")) +
  geom_smooth(aes(x=i,y=real_i*1e3,color="in well"),method = lm,size=.5,se=F) +
  geom_smooth(aes(x=o,y=real_o*1e3,color="out well"),method = lm,size=.5,se=F) +
  theme_tsj() +
  labs(x = "Dark Corrected DN/ms",
       y = bquote("Radiance [mW"~m^{-2}~nm^{-1}~sr^{-1}~"]"),
       color="") +
  geom_text(aes(x=1500,
                y=10,
                label = paste0(
                  "y = ",
                  as.character(round( in_well_a*1e3, 4)),
                  "x + ",
                  as.character(round( in_well_b*1e3, 3))),
                color="in well")) +
  geom_text(aes(x=1500,
                y=20,
                label = paste0(
                  "y = ",
                  as.character(round( out_well_a*1e3, 4)),
                  "x + ",
                  as.character(round( out_well_b*1e3, 3))),
                color="out well"))

fit <- lm( real_i ~ i, df)
in_well_b <- fit$coefficients[1]
in_well_a <- fit$coefficients[2]  

fit <- lm( real_o ~ o, df)
out_well_b <- fit$coefficients[1]
out_well_a <- fit$coefficients[2]

  
#cosine correction:-----------------------------------------------
#cosine_df <- read.csv("cal_files/QEP02723_cosine_calibration_20200428-0124.csv")
#cosine_df$pan_angle <- cosine_df$pan_angle - cosine_df$pan_angle[ cosine_df$DN_in == max( cosine_df$DN_in) ]
#cosine_df$corr_DN_in <- cosine_df$DN_in/max(cosine_df$DN_in)
#cosine_df$corr_DN_out<- cosine_df$DN_out/max(cosine_df$DN_out)
#ggplot( cosine_df, aes(x=pan_angle) ) +
#  geom_point( aes(y=corr_DN_in, color="well") ) +
#  geom_line( aes(y=corr_DN_in, color="well") ) +
#  geom_point( aes(y=corr_DN_out, color="shoulder") ) +
#  geom_line( aes(y=corr_DN_out, color="shoulder") ) +
#  theme_tsj() +
#  labs(x="Pan Angle [degrees]", y= "Fraction of Signal",color="Pixel" )



#sky_df$i <- sky_df$i/approx( cosine_df$pan_angle, cosine_df$corr_DN_in, sky_df$SZA)
#sky_df$o <- sky_df$o/approx( cosine_df$pan_angle, cosine_df$corr_DN_out, sky_df$SZA)

df$i_sky <- approx( sky_df$time, sky_df$i, df$time )$y
df$o_sky <- approx( sky_df$time, sky_df$o, df$time )$y
#df$i_spect <- approx( spect_df$time, spect_df$i , df$time )$y
#df$o_spect <- approx( spect_df$time, spect_df$o , df$time )$y

#calculate SIF:
df$SIF_sky   <- ( df$i*df$o_sky - df$o*df$i_sky) / (df$o_sky - df$i_sky)
#df$SIF_spect <- ( df$i*df$o_spect - df$o*df$i_spect) / (df$o_spect - df$i_spect)

plt_all_sif <- ggplot(df, aes(x=hour,y=SIF_sky,color=target))+
  geom_point(size=.4) +
  facet_wrap(~date,nrow=2) +
  ylim(0,.3)

#sdf <- df[ df$target %in% c("sky","spect","grass_1","grass_2") , ]
#sdf <- df[ df$okta < 2, ]
df$sample_no <- 0
for( target in targets){
  sdf <- df[ df$target == target , ]
  df$sample_no[ df$target == target] <- 1 + c(0, cumsum( diff(sdf$time) > 70 ) )
}



#sdf <- df[ (df$time > as.POSIXct("2020/05/07 19:00:20", tz="UTC") ) & 
#           (df$time < as.POSIXct("2020/05/07 19:02:40", tz="UTC") ) , ]


#sdf$time <- round_date( df$time, unit="10 mins")
noise_df <- aggregate( .~ sample_no+target+date,df,"sd" )
sdf <- aggregate( .~ sample_no+target+date,df,"median" )
n_df <- aggregate(.~ sample_no+target+date,df,"length")
sdf$SIF_sky_sd <- noise_df$SIF_sky
sdf$n_samples <- n_df$SIF_sky
#sdf$SIF_spect_sd <- noise_df$SIF_spect

#sdf <- sdf[ sdf$SIF_sky_sd < .001 , ]

#sdf <- df[ df$sample_no == 205, ]

bar_df <- aggregate( .~target+date,sdf,"mean")
sdf2 <- sdf[ sdf$sample_no ==1, ]
ggplot( bar_df,aes(x=pan,y=tilt,color=SIF_sky)) +
  geom_point(size=6) +
  xlim(-30,-12) +
  ylim(-25,-10) +
  theme_tsj() +
  scale_color_gradientn(colors=lumen) +
  labs(title=as.POSIXct(sdf$time[1],origin="1970-01-01") )
  

pal <- c( '#44AA88', #Grass
          '#668822', #maple
          '#B9ABA7',  #pavement
          '#5566AA', #river
          '#555555', #roof
          '#9D704A',#tilia
          '#668822',
          '#99BB55')  #tree

sdf$target[ sdf$target == "pavememt"] <- "pavement"
sdf$target[ sdf$target == "grass_1"] <- "grass"
sdf$target[ sdf$target == "maple_1"] <- "maple"
sdf$target[ sdf$target == "tilia_1"] <- "tilia"

sdf <- sdf[ sdf$target != "sky", ]

ggplot( sdf, aes(x=hour,y=SIF_sky,color=target)) +
  geom_point(alpha=.5) +
  geom_smooth(method = "loess", se=F) +
  scale_color_manual(values=pal) +
  ylim(-.001,.025) +
  #facet_wrap(~date,ncol=1) +
  scale_x_continuous(breaks=seq(12,24,.5),limits=c(14,15.5)) +
  theme_tsj() +
  labs(x="Hour of Day [UTC]", y=bquote("SIF [mW"~m^{-2}~nm^{-1}~sr^{-1}~"]"), color="Target")


ggplot( df, aes(x=hour,y=SIF_sky)) +
  geom_line() +
  #ylim(-.01,.08) +
  facet_wrap(~target) +
  scale_x_continuous(breaks=seq(12,24,2),limits=c(10,24)) +
  theme_tsj() +
  labs(x="Hour of Day [UTC]", y="SIF", color="Target")






